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We present a spectral method for solving elliptic equations which arise 
in general relativity, namely three-dimensional scalar Poisson equations, as 
well as generalized vectorial Poisson equations of the type AN + AV(V ■ 
N) = S with A / -1. The source can extend in all the Euclidean space R 3 , 
provided it decays at least as r -3 . A multi-domain approach is used, along 
with spherical coordinates (r, 9, <f>). In each domain, Chebyshev polynomi- 
als (in r or 1/r) and spherical harmonics (in 9 and <j>) expansions are used. 
If the source decays as r~ k the error of the numerical solution is shown 
to decrease at least as N~ 2(k ~ 2 \ where N is the number of Chebyshev 
coefficients. The error is even evanescent, i.e. decreases as exp(— N), if the 
source does not contain any spherical harmonics of index / > k — 3 (scalar 
case) or I > k — 5 (vectorial case). 
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1. INTRODUCTION 
1.1. Scalar and vectorial Poisson equations with non-compact sources 

The most common elliptic equations which occur in numerical relativity (for a 
recent review see [1]) are the scalar Poisson equation 

= S , (1) 

and the (generalized) vector Poisson equation 

AN + AV(V ■ N) = S , (2) 

where A is a constant different from —1, typically A = 1/3. Contrary to the Newto- 
nian case, where the source term S contains only the matter density, the sources of 
these equations have a non-compact support. Moreover, the Einstein equations be- 
ing non linear, the sources S and S depends (usually quadratically) on the solutions 
F and N. This means that equations (1) and (2) must be solved by iterations. 
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More precisely, within the 3+1 formalism (also called Cauchy formulation) of gen- 
eral relativity (see [2] for a review), the 10 Einstein equations can be decomposed 
into a set of 6 second order evolution equations and 4 constraint equations: a scalar 
one, the so-called Hamiltonian constraint, and a vectorial one, the so-called mo- 
mentum constraint (see [3] for an extensive discussion of the constraints equations). 
The PDE type (i.e. hyperbolic, parabolic or elliptic) of these equations depend 
on the coordinates chosen to describe the space-time manifold. Let us recall that 
within the 3+1 formalism, the space-time is foliated in a family of space- like slices 
S t , labeled by the time coordinate t. The space-time 4-metric is then entirely de- 
scribed by the induced 3-metric 7^ of the hypersurfaces £4 along with the extrinsic 
curvature tensor Kij of S t . 

In this context, a typical example for Eq. (1) is the equation for the lapse function 
for the choice of time coordinate corresponding to a maximal slicing of space-time 1 
(sec e.g. [4]). Another example is provided by York treatment of the initial- value 
problem of general relativity [5], according to which the Hamiltonian constraint 
equation results in an elliptic equation of the type (1) for the conformal factor of 
the spatial metric 7^, with a term F~ 7 in S. 

Regarding the vector Poisson equation (2), it also appears in York formulation of 
the initial-value problem for the vector which enters in the longitudinal part of the 
transverse-traceless decomposition of the extrinsic curvature tensor K^. Indeed 
the momentum constraint determines the longitudinal part of according to the 
equation 2 



where Vj is the covariant derivative associated with the 3-metric 7^ , J 1 is the mat- 
ter momentum density, and maximal slicing is assumed [K\ — 0). More generally, 
the vector Poisson equation (2) with A = 1/3 occurs each time one has to perform 
the transverse-traceless decomposition of a symmetric tensor field T y defined on 
a Riemannian three-manifold with metric 7^ . Following [5, 6], this decomposition 
writes 



VjK 13 = 8wJ l , 



(3) 




where T = "fuT kl . T^ T is the transverse-traceless part, (LY) 13 the longitudinal 
trace-free one and |T7 y the trace part. The longitudinal part is expressible in 
term of a vector Y l , by means of the conformal Killing operator : 




Performing the decomposition reduces to the finding of the vector field Y. Consid- 
ering the divergence of Eq. (4) , Y appears to be the solution of the equation 




1 This Poisson equation for the lapse function reduces to the usual Poisson equation for the 
gravitational potential at the Newtonian limit 

2 Einstein convention of summation on repeated indices is used 
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where Rj is the Ricci tensor associated with the metric 7^ . This is a vectorial 
Poisson equation of type (2) with A = 5 (involving the so-called conformal Laplace 
operator). Let us mention that, in the general case, it must be solved by iteration 
for Y is present in the source term. 

Another example for the vectorial Poisson equation (2) is provided by the so- 
called minimal distortion [4] choice of coordinates in the spatial hypersurfaces E t . 
The unknown vector TV is in this case the shift vector which defines the propaga- 
tion of the spatial coordinates x % from one slice S t to the next one T, t +dt- It is 
this vectorial Poisson equation, which is a special form of Eq. (2) with A = |, that 
originally motivated our study of this subject. Let us mention that the confor- 
mal Killing operator and the associated vectorial Poisson equation also appear in 
the "thin-sandwich" formulation, where the spatial geometry is given on two close 
hypersurfaces (see [3, 7] for more details). 

1.2. Treatment by means of spectral methods 

Solving elliptic equations is often considered as a CPU-time consuming task. 
Spectral methods [8, 9] seems attractive in this respect because they provide accu- 
rate results with reasonable sampling, as compared with finite difference methods 
for example. We refer the interested reader to Rcfs. [10, 11] for a review of the use of 
spectral methods in relativistic astrophysics. Let us simply mention here that our 
group has previously developed a spectral method, using Chcbyshcv polynomials 
and spherical harmonics to solve three-dimensional scalar Poisson equations with a 
compact source [12]. However as recalled above, the elliptic equations which arise 
from numerical relativity have non-compact sources. This means in particular that 
the infinity is the only location to impose exact boundary conditions (flat space- 
time). In order to tackle this, we have introduced a multi-domain approach [13] 
within which the last domain extends up to infinity, thanks to some compactifica- 
tion. This approach has another nice feature, for it is avoiding Gibbs phenomena: 
a physical discontinuity can be located at the boundary between two domains, so 
that all the considered fields are smooth in each domains. 

In this article, we extend the single-domain spectral method for the scalar Poisson 
equation (1) presented in [12] to the multi-domain case, which enables in particular 
to treat non-compact sources provided they decay at least as r~ 3 when r — > 00. 
Based on this scalar Poisson solver, we treat the generalized vectorial Poisson equa- 
tion (2). We consider three different schemes proposed in the literature to reduce 
the resolution of (2) to 4 scalar Poisson equations, namely the schemes of Bowen & 
York [14], Oohara & Nakamura [15] and Oohara, Nakamura & Shibata [16]. These 
schemes have been originally implemented on finite (single) domains and with fi- 
nite difference methods. We study here their applicability to infinite domains and 
spectral methods. 

The solvers presented in this work deal with three-dimensional flat spaces where 
V denotes the ordinary derivation. More general cases (i.e. Laplacian operator 
associated with a curved metric) can be solved by iteration. In all the following 
we will assume that there exists a unique solution of both the scalar and vectorial 
equation that is C°° by parts, C 1 everywhere and that is going to zero at infinity. For 
known results about the existence and uniqueness of solution of partial derivative 
systems see for example [17]. 
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This paper is organized as follows. In Sec. 2 we present the numerical scheme 
used to solve the scalar Poisson equation with our multi-domain spectral method. 
This scheme is tested in Sec. 3 using comparison with analytical solutions of var- 
ious behaviors. This study leads us to establish the convergence properties of the 
algorithm. Sec. 4 is devoted to the study the three different schemes mentioned 
above to solve the vectorial Poisson equation (2). As for the scalar Poisson equa- 
tion, the implemented schemes are tested in Sec. 5 and their convergence properties 
exhibited. In Sec. 6 we give to the reader some indications about some extensions 
of this work that have been successfully conducted or under work. 

2. SCALAR POISSON EQUATION 
2.1. Spectral expansions 

As described in previous articles [10, 12], spherical coordinates (r, 9, <j)) are used; 
the fields are expanded in spherical harmonics Y™ (9, </>) and a Chebyshev expansion 
is performed with respect to the r coordinate. Doing so the resolution of the scalar 
Poisson equation is reduced to find, for each couple (l,m), the solution of 

where / and s are functions of r solely, being respectively the coefficients of YJ m in 
the solution F and in the source S. 

f and s are expanded in Chebyshev polynomials (hereafter referred as for the 
polynomial of order i) so that the inversion of the operator on the left-hand-side of 
Eq. (7) is reduced to a matrix inversion. 

As recalled above, the present work improves that presented in [12] for we are 
allowing a source that is not compactly supported. To take care of this, we will 
divide space in three type of domains, following [13] 

• One kernel, a sphere centered at the origin and being the only domain consid- 
ered in [12]. In such a domain r is given by r — ax, where x G [0, 1], with a > 0. 
The functions are expanded in Chebyshev polynomials in x with a definite parity 
to ensure regularity at the origin : only even (rcsp. odd) polynomials arc involved 
for I even (resp. odd). 

• An arbitrary number, including zero, of shells, domains where r = ax + (3, 
x G [—1,1]. We have the following conditions : a > and j3 > a, so that r is 
increasing with x and never equal to zero. In the shells, the functions are expanded 
in usual Chebyshev polynomials, with no parity requirement. 

• One external domain, extending to infinity, where r is given by u = r _1 = 
a(x — 1), a being negative and x G [—1, 1]. Once more the functions are given as a 
sum of Chebyshev polynomials in x. 

2.2. The matrices 

Before doing any operator-inversion, one has to take care about singularities at 
the origin and at infinity. For example, because of division by r 2 , the solution of 
the equation, must be decreasing as r 2 at the origin to be associated with a non- 
singular source. We choose to treat that by subtracting finite parts of the solution 
at the point of singularity. 
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Before describing that more precisely, let us mention another method for solving 
that problem, presented in [12]. In this Reference, the functions are expanded 
on a new set of basis- functions, that verify individually the regularity conditions 
(Galerkin basis) . For example, + Tj is used in the kernel, making all the 
basis-functions decrease as r 2 at the origin. 

• In the kernel, we have to take care of a singularity at the origin due to the 
division by r 2 . To avoid this we construct an operator without the finite part of / 
at x = 0. Thus the operator is, expressed in terms of x, 

the source s being multiplied by a 2 . 

r 2 

• In the shells there is no singularity, so we can multiply the source by — ^ and 

a z 

invert the following operator 



• In the external domain Eq. (7), once re- written in terms of u = -, becomes 



i4 /d 2 / _ 1(1 + 1) \ 

\du 2 u 2 ' ) 



(10) 



We consider the three following possibilities. 



— First multiplying the source by r 4 in the external domain, a singularity 
occurs at r = oo, that is x = 1. We handle it like in the kernel, by subtraction of 
the finite part of / in 1 , and we use the following operator 

— If the source is multiplied by r 3 , a singularity occurs r — oo, that is x = 1 
and is handled by the finite part method, so that the operator becomes 

^-"-"ff-Flf-'f 1 "- (12) 

— If the source is only multiplied by r 2 , we invert the non-singular operator 

A/ = Or-l) 2 0-Z(Z + l)/. (13) 

In all cases we have to multiply s by a 2 . Let us stress out that those three operators 
are not fully equivalent in actual physical calculations based on iterative schemes. 
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The effective source (i.e. r k S) being given, the solution will have less high-frequency 
terms (Chebyshev polynomials of high-order), if the number k is high. Those high- 
frequency terms may cause instabilities in an iterative procedure, so we always use 
the r 4 S scheme except for a source decreasing like r 3 S at infinity. 



As an illustration, here is the matrix constructed in the kernel, with I — 2 and 9 



r (Chebyshev polynomials To, 


T 2 , 


. Tie) 




(0 





56 


96 


304 


480 


936 


1344 


2144 








56 


240 


472 


1056 


1656 


2832 


3992 











144 


432 


848 


1632 


2512 


3984 














264 


688 


1320 


2336 


3528 

















416 


1008 


1888 


3168 




















600 


1392 


2552 























816 


1840 


























1064 


\o 



























2.3. The banded-matrices 

The constructed matrices are not suitable for numerical purpose. The inversion 
would be much more rapid and much more efficient if we could work on band- 
matrix instead of triangular ones. The operators being second order operators on 
a set of orthogonal functions, there must exist a linear combination of the lines so 
that the matrices are reduced to banded-ones (see [18]). 
We exhibit here the combination we used in each domain : 

• In the kernel, the Chebyshev polynomials are either odd or even, depending on 
the parity of I. The combination is independent of the actual value of I except for 
its parity. 

When the Chebyshev polynomials are even we use : 

U = (1 + 51) U - L t+2 (for < i < N - 3) (14) 
U = Li- L l+2 (for < i < N - 5) (15) 
Li = U - L l+1 (for < i < N - 5) (16) 



and when they are odd : 



U 


= Li 


— Li+2 


(for 


< i < N 


-3) 


(17) 


Li 


= Li 


— Li + 2 


(for 


< i < N 


-5) 


(18) 


Li 


= U 


— Li+i 


(for 


< i < N 


-5) 


(19) 



where Li denotes the line number i and N is the number of Chebyshev polynomials 
involved in the expansion. 

In both cases the resulting matrix is a 4-band one. 
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• In the shells, the basis of decomposition contains all the Chebyshev polynomi- 
als. The combination is the following : 



(l + 5q) Li - L i+2 



i + 1 



(for < i < N - 3) 



Li = Li - L i+2 (for < i < N - 5) 



(20) 
(21) 



The resulting matrix is a 5-band one. 

• In the external domain, the combination depends of the type of the constructed 
operator. 

— If the source is multiplied by r 4 the combination is : 



Li 


= (1 


+ aj) Li 


- L i+2 (for < i < N - 3) 


(22) 


Li 


= Lt 


— Li+2 


(for < i < N - 5) 


(23) 


L'i 


= U 


— Li+i 


(for < i < N - 5) 


(24) 


Li 




~ L'i+2 


(for < i < N - 5) 


(25) 



The resulting matrix is a 4-band one. 

— If the source is multiplied by r 3 the combination is : 

Li = (l + S z ) Li - L t+2 (for < i < N 
U = Li - L t+2 (for < i < N - 5) 
Li = Li + L l+l (for < i < N - 5) 



3) 



(26) 
(27) 
(28) 



The resulting matrix is a 4-band one. 

— If the source is only multiplied by r 2 , the combination is the same than the 
one used in the kernel for even polynomials. Then, the resulting matrix is a 6-band 
one. 

Of course to maintain the solution, the same linear combination is performed on 
the coefficients of s. 

The banded matrix associated with the one presented above (in the kernel with 
I = 2 and N = 9) is : 

\ 








56 


-336 


-200 




















56 


96 


-488 


336 




















144 


168 


-672 


-504 




















264 


272 


-888 


-704 




















416 


408 


-1136 


-2000 




















600 


1392 


1488 























816 


1840 


























1064 


\o 



























2.4. Homogeneous solutions 

Due to the presence of homogeneous solutions, the banded-matrices are not in- 
vertible. The operator given by Eq. (7) has two homogeneous solutions which are 
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r l and r~( l+1 \ Those functions are eigenvectors of the matrix with the eigenvalue 
0. In the kernel and the external domain, the use of the finite part of the solution 
can sometimes introduce other homogeneous solutions. 
Let us summarize the number of such eigenvectors in each case : 

• In the kernel the solution in r~(' +1 ) is singular for r = and so is not taken 
into account. We have one additional homogeneous solution, arising from the finite 
part : T for I even and T\ for I odd. 

The parity of the Chebyshev polynomials is the same than that of I so the eigen- 
vectors are : 

- T only for I = 0. 

- Ti only for I = 1. 

- r l and T for I > 2, even. 

- r l and 7\ for I > 3, odd. 

• In the shells we have to take into account the two usual homogeneous solutions, 
which are not singular in this case. We could remark that if r l is exactly described 
by the Chebyshev expansion, and so implies an exact zero determinant for the 
matrix, it is not the case for the fractional solution r~^ l+1 \ This one is not given 
by a finite sum of Chebyshev polynomials but rather by an infinite sum implying 
the result would be worse and worse as the number of coefficients increases, for the 
determinant of the matrix would be closer and closer to 0. So, to deal with this, 
we have to take into account that the eigenvalue is of order 2 even if r~(' +1 ) only 
becomes an exact eigenvector for an infinite number of coefficients. 

• In the external domain, the solution r l is singular at infinity except for I = 0. 
r -(l+i) j g a i wa y S acceptable. 

If the source is multiplied by r 4 , the finite part introduce two others eigenvectors 
of eigenvalue : T and T\. So the situation is : 

- T and T x for I = 0. 

- T , Ti and r - {l+ ^ for I > 1. 

If the source is multiplied by r 3 , the finite part only introduce one other eigenvector 
of eigenvalue : To, and the situation is : 

- T and r-(' +1 ) for all I. 

If the source is multiplied by r 2 , they are no other solutions than the usual ones 
which gives : 

- T and T x for I = 0. 

- r-(' +1 ) for I > 1. 

From the above discussion we are able to determine the order p of the eigenvalue 
0. The banded-matrices are then amputated from their p first columns and their p 
last lines resulting in invertible banded-matrices. We abandon the p last coefficients 
of the source. Doing so, we find a particular solution of the system which has its p 
first coefficients undefined and thereafter set to zero. 
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In particular the previously presented matrix (in the kernel, for Z = 2 and N = 9) 
becomes : 



/56 


-336 


-200 











° \ 


56 


96 


-488 


336 














144 


168 


-672 


-504 














264 


272 


-888 


-704 














416 


408 


-1136 


-2000 














600 


1392 


1488 


Vo 














816 


1840 / 



Before solving the system, an LU decomposition is performed using LAPACK 
(Linear Algebra PACKage) [19] for purpose of rapidity. LAPACK is also used for 
the resolution of the system. 

2.5. Regularity and boundary conditions 

In this section we will show how some homogeneous solutions are used to main- 
tain regularity and satisfy the boundary conditions. We will concentrate on the 
boundary condition / = at infinity. 

• In the kernel, the operator is singular only for I > 2. If it is the case, to 
maintain regularity, / has to verify the following conditions : 

/(0) = (29) 
f(0) = (30) 

Thanks to the parity of the Chebyshev expansion, one of these conditions is always 
fulfilled, depending on the parity of So we perform a linear combination of the 
solution with either T or 7\ to fulfill the other one. Nothing has to be done for 
I < 1. 

• In the shells, nothing has to be done for there arc neither boundary conditions 
nor singularities. 

• In the external domain we should, once more, discriminate between three cases 

— If the source is multiplied by r 4 , we must impose / (1) = to satisfy the 
boundary condition ; this is done by performing a linear combination of the solution 
and T . Then for I > 1, for reasons of regularity, we must have /' (1) = 0, a 
condition which is obtained by linear combination with T\. 

— In the case of a source multiplied by r 3 , the boundary condition f (1) = 
ensures regularity ; we impose it by performing a linear combination of the solution 
and T . 

— If the source is multiplied by r 2 , the situation is a bit more subtle. There 
is no condition of regularity, but the boundary condition imposes that / (1) = 0. 
Then one can show that this implies that the source decreases as r -3 at infinity. 
Conversely, if the source decreases as r~ 3 , it implies, for / ^ 0, that / (1) = 0. 

So to verify boundary conditions we only consider sources decreasing as r~ 3 . It 
implies that the boundary condition is automatically verified for I 0. We only 
impose it for I = 0, by doing a linear combination of the solution and Tq. 
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Let us emphasize that this is only the theoretical aspect of the problem. During 
an actual physical calculation, the sources of the Poisson equation are themselves 
numerically given so that they might, due to computational errors, not decrease 
exactly like r~ 3 . In such a case one should be cautious, for the solution will not 
exactly be zero at infinity. A possible treatment is to enforce the r~ 3 decay by 
slightly modifying the source prior to the resolution of the Poisson equation. 

2.6. Continuity 

At this stage, for each (I, m), we are left with a particular solution in each domain, 
one homogeneous solution in the kernel and in the external domain, and two in each 
shell. The last linear combinations will be performed to ensure the continuity of 
the solution and of its first derivative across each boundary. 

The simplest case is when the angular sampling is the same in every domain 
(i.e. the same numbers of point in 9 and <fi). The unknowns are the coefficients of 
the homogeneous solutions in the physical solution and the equations are given by 
matching / and its derivative across each boundary. It is easy to see that there 
is exactly the same number of equations and of unknown quantities, resulting in a 
uniquely determined solution. 

If the angular sampling is not the same, the situation is a bit more complex, for 
some Y; m may not be present in some domains. At each boundary, for each (/, m), 
three situations can occur : 

• the harmonic is present in the two domains : we perform the matching of both 
/ and its derivative. 

• the harmonic is present in neither domains : no equation is written. 

• the harmonic is present only in one domain : we assure the continuity of / 
supposing that the harmonic has its coefficient equal to in the domain where it 
is not present. We perform no matching for its derivative. 

This procedure results in a system of equations that admit a unique set of solutions. 
We have imposed exactly as much continuity as the sampling allowed us. 
To illustrate this, let us take the situation given by Tab. 1, for a specific value 
(Z,m) : 

TABLE 1 



Example of the situation before doing the connection across each boundary 



Domain 


Bounds 


yrn 


Particular 
solutions 


Homogeneous 
solutions 


Unknowns 





< r < Ri 


Yes 


fo 


r l 


Cto 


1 


Ri < r < R 2 


Yes 


fi 


r l and r^ 1 ' 


ai and /3i 


2 


R 2 <r <R 3 


No 








3 


R 3 <r <R 4 


No 








4 


Ri < r < oo 


Yes 


u 


r -G+i) 


134 



In that situation the domain is the kernel and the domain 4 is the external 
compactificd region. The column labeled Y" 1 denotes the presence or the absence of 
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the considered spherical harmonic in each domain. The particular and homogeneous 
solutions are expressed taking into account the sampling and the nature of each 
domain. The unknowns are the coefficients of the homogeneous solutions labeled 
a for r l and (3 for r~( L+1 \ Using the procedure described above we obtain the 
following equations : 

• For r — Ri, the spherical harmonic is present in both domains, so we have to 
write the continuity of the solution and its derivative, which gives 

f (Ri)+a R[ = fi(Ri)+aiR[+PiRi {l+1) (31) 
f^(Ri)+laoR ( t 1) = f[(Ri) + laiR ( t 1] -(l + l)PiRi {l+2) . (32) 

• For r = R 2 , the spherical harmonic is present only in the domain 1 and so we 
write only the continuity of the solution assuming that it is zero in the domain 2 

fi (R 2 ) + aiR l 2 + [3iR- (l+1) = 0. (33) 

• For r = i?3, no equation is written, for the harmonic is absent on both sides of 
the boundary. 

• For r = Ri, the situation is the same as at r = R 2 



U(R4) + faR? l+L> = 0. (34) 

We have now four independent equations which are solved to find the unknowns 
oeo, oil, Pi and (3^. 

All this procedure enables us to find a unique solution of the scalar Poisson 
equation, solution to be regular everywhere, continuous, like its derivative, and 
that is zero at infinity. We should point, once more, that the source must decrease 
at least as r -3 , for this to be possible. 

3. CONVERGENCE PROPERTIES OF THE SCALAR POISSON 

EQUATION SOLVER 
3.1. Position of the problem 

We wish to study the convergence of our algorithm, depending on the number 
of coefficients chosen for the r-expansion. The number of points for and <j> does 
not change the precision of the result, as long as we have enough points, that is 
enough spherical harmonics, to describe the source properly. However, concerning 
r, we perform matrix-inversion and we expect a better precision as the number of 
coefficients increases. 

It is well known (see [8, 9]) that with spectral method, the error is evanescent, 
i.e. decreasing as exp(— N), N being the number of coefficients, as long as we 
are working with functions that are C°°. If the functions are only C p , the error is 
decreasing as _/V~( p+1 ) solely. This is known as the Gibbs phenomenon. 

The various domains of our multi-domains method [13] are intended to fit the 
surfaces of discontinuity, for example the surface of a star (see [20] for an application 
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to stars with discontinuous density profiles, like strange stars). Doing so, each 
function is C°° in each domain, removing any Gibbs phenomenon. 

To test the validity of our numerical scheme, we compare calculated solutions to 
analytical ones. We estimate the relative error as the infinite norm of the difference 
over the infinite norm of the analytical solution. We will present some examples for 
the construction of analytical solutions. Using only C°° functions we expect errors 
to be evanescent. 

But this is not so simple. It can be easily shown that, generally, the particular 
solutions obtained by the inversion of the operator for each (/, m) are of polynomial 
or fractional type. Those functions are exactly described by Chebyshev polynomials 
inrorr -1 . This is true except in two cases, related to the homogeneous solutions : 

• a source in r l ~ 2 , will give rise to a particular solution in r l lnr. 

• a source in r~" +3 > , will be associated with a particular solution in r~(' +1 ) lnr. 

In such cases, we expect some problems, for the description of logarithm functions 
in terms of Chebyshev polynomials may not be accurate. To be more precise about 
this effect, let us study the situation in each type of domain. 

3.1.1. In the kernel 

In the kernel and for reason of regularity, sources in r~(' +3 ) are obviously never 
present. At first sight, the case of a source in r l ~ 2 seems to be more problem- 
atic ; but let us recall that this source has to be the factor of Y™ . Can we have 
a source containing terms like Y^r 1-2 ? To answer this question we refer to [10] 
where it is shown that, for a regular function (i.e. function expandable as a poly- 
nomial series in Cartesian coordinates (a;, y, z) associated with (r, 9, </>)), terms like 
r aym are p resen ^ m spectral expansion only if a > I. So, sources leading to a 
In function in the kernel are not regular at the origin. To conclude, we expect no 
problem connected with particular solutions containing In functions in the kernel, 
at least with physical regular sources. However, let us mention the fact that if the 
source is the result of some calculation, it might contain some unphysical terms 
due to computational errors. Those terms might give rise to some logarithmic 
functions. 

3.1.2. In the shells 

As usual there are no regularity prescriptions in the shells. The two types of 
particular solutions can appear. To investigate more precisely the effects of the 
logarithm, we studied the behavior of the error performed by expanding the two 
types of particular solutions in Chebyshev polynomials. 

We constructed the two following exact particular solution r l lnr and r~(' +1 ) lnr, 
approached them by a sum of Chebyshev polynomials in x, the relation between r 
and x being r — ax + (3. Then, we estimate the error by the same method as the 
one described before. 

Fig. 1 shows an evanescent error. The functions containing logarithm are thus 
rather well described in a shell. This is due to the fact that the In functions are 
bounded in such domains and not going to infinite values. More precisely, the 
r'lnr and r~^ +1 ^lnr functions are C°° in the shells so that the error should be 
evanescent. Let us mention that this result does not depend on the choice made 
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FIG. 1. Relative difference (infinite norm) between the particular solutions with logarithm 
and their truncated Chcbyshev expansions, in a shell. 
The scale for the number of coefficients is linear. 

The solid lines represent the r l lnr functions and the dashed lines the r~ lnr ones. 
The circles represent the case 1 = 0, the squares 1 = 1 and the diamonds 1 = 2. 
This plot has been obtained using a = 0.5 and (3 = 1.5. 



for a and (3. To conclude we expect no problem to rise from the presence of such 
particular solutions in the shells. 

3.1.3. In the external compactified domain 

The particular solution in r l In r is not going to zero at infinity, and so cannot 
appear in the external domain. But the other type of particular solutions r~^ l+v> In r 
is likely to appear. We investigate their effect by the same method than the one used 
in the shells, that is determining the behavior of the error done by interpolating 
the exact solution by a finite sum of Chebyshev polynomials. 

Fig. 2 shows that the error is no longer evanescent but follows a power law. 
The error is decreasing faster and faster as I increases, for the associated particular 
solution is being better approached by Chebyshev polynomials. In other words, the 
function r~(' +1 ) lnr = — u^ +1 ^ lnu is not C°°, for its (I + l) th derivative contains 
terms in i, not regular at spatial infinity, that is for u = 0. More precisely, Fig. 3 
shows the value of the exponent as a function of I. 

We can conclude that the error done by expanding r~^ +1 ^lnr in Chebyshev 
polynomials follows a power-law and that it is decreasing faster than N~ 2 ( l+1 \ Wc 
will use this to explain some features of our scalar and vectorial Poisson equation 
solvers. 



3.2. Accuracy estimated by comparison with analytical solutions 

From the results of the previous section, we expect an evanescent error for the 
resolution of the scalar Poisson equation when there is no particular solution con- 
taining any logarithm in the external domain and an error following a power-law 
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FIG. 2. Relative difference (infinite norm) between the particular solutions with logarithm 
and their truncated Chebyshev expansions, in the external domain. 
The scale for the number of coefficients is logarithmic. 

The circles represent the case I = 0, the squares I = 1 and the diamonds / = 2. 
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FIG. 3. Exponent of the power-law followed by the error shown in Fig. 2, as a function of 

I. 
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FIG. 4. Error on the resolution of the scalar Poisson equation for a spherically symmetric 
source extending to infinity. 

The solid lines represent the r 4 S scheme, the dotted lines the r 3 S scheme and the dashed lines 
the r 2 S one. 

The scale for the number of coefficients is linear. 

The circles represent the error in the kernel, the squares in the shell and the diamonds in the 
external domain. 

when such solutions appear. We present here some results that illustrate this be- 
havior and lead to two properties about the error. 

3.2.1. Spherically symmetric source 

First of all, let us consider a simple case for which we do not expect any Gibbs- 
like phenomenon : a spherically symmetric source decreasing as r~ 4 . In fact, the 
only harmonic present in this source is I = 0, that would imply a In solution only 
for a source in r~ 3 . We choose a source S decreasing as r~ 4 in the external domain 
and a polynomial one, such that the solution is not singular in the kernel. The 
associated solution F can be found analytically. 
In the external domain, for r > R, we have 

s = - ; F = ^-^V (35) 



r 

and for r < R, 



S = R-- ■ F= — - — --R i . (36) 
R ' 6 20R 4 v ' 

As expected Fig. 4 shows an evanescent error, with some saturation at the level 
of 10 -15 due to the round-off error, the calculation being performed in double 
precision. No significant difference can be seen between the three schemes. 



3.2.2. Compact source 
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FIG. 5. Error on the resolution of the Poisson like equation for a non-spherical compact 
source with 1 = 2. 

The scale for the number of coefficients is linear. 

The circles represent the error in the kernel, the squares in the shell and the diamonds in the 
external domain. 

Another interesting case is that of a source with a compact support, that is a 
source which is zero in the external domain. As for the previous case we do not 
expect any Gibbs phenomenon. In the external domain let us choose the following 
analytical solution 



F = Y, 



1 ' 



(37) 



This solution leads to a source that vanishes in the external domain. To avoid any 
singularity at the center, we choose the later function as a solution of the equation 
for r < R 



F = Y, 



(21 + 5) 



2jR 2/+3 



(21 + 3) 



2 .R2i+5 



(38) 



This solution has been chosen so that F and its first derivative with respect to 
r are continuous at r — R, properties of the solution given by our algorithm. The 
associated source, for r < R, is found by taking the Laplacian of F 



(21 + 5) (21 + 3) 



1 



R 2l+3 



(21 + 3) (41 + 10) 



R 2l+5 



(39) 



So we constructed a non-spherically symmetric compact source, which contains 
only one spherical harmonic. We chose for simplicity m — 0, for we do not expect 
any variation with to, the later being absent of the inverted operator. 

As expected, Fig. 5 shows an evanescent error down to a saturation value of 
approximatively 10~ 14 . 
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3.2.3. A logarithm in a shell 

The last case with an evanescent error we considered is the one where the prob- 
lematic particular solutions (i.e. containing a logarithm) appear only in a shell 
bounded by R\ < r < R2. We choose a source s that implies the appearance of 
both type of particular solutions. Let F be the associated solution. In the shell, 
for Ri < r < i?2, we have 

1 z 2 

Sshcll = -^+3^-1 (4 
lnr ln.Ri-1 1 ( z 2 \/l 2l /hii? 2 1\ 2 #1 1 

F ^ = -- + ^— + r 2 + { 3 ^- 1 ) [r l * r -{— + 25) r + 25^ 

For simplicity, we take S = in the kernel and in the external domain, the 
solution being chosen, once more, by continuity across the boundaries 

lnfli -lnR 2 1 R\-Rl ( 7?_ \ 
r r 3 25 \ r 2 J 



-^kernel 
-^external — 



The result presented in Fig. 6 shows an evanescent error, confirming that the 
presence of logarithm function is only a problem in the external domain. Once 
more let us mention that this is due to the fact that the logarithm functions are 
bounded in a shell and not going to infinite values. Such bounded functions are 
rather well described in terms of Chebyshev polynomials. 

3.2.4- The Gibbs phenomenon 

Let us now consider a case where the particular solution contains a logarithm in 
the external domain. Following the construction of the source and solution of Sec. 
3.2.2 let us take the following source in the external domain 



S=-Yt°^, (42) 



and £ = for r < R. The associated unique solution is 



J 



F = Y,° 5 for r<R (43) 

(21 + l) 2 R 2l+1 

n ln(r) - In (i?) + ^yi-r 

F = Y?—^— \ '1 21+1 for r > R. 

1 [21 + 1 )r l + 1 



Fig. 7 presents an example of the obtained result for each of the three schemes 
discussed in Sec. 2.2. A logarithm being present in the solution, the error is no 
longer evanescent and follows a power-law. One important feature is that the r 4 S 
scheme is converging much less rapidly than the r 3 S and r 2 S ones. It may come 
from the fact that for a given source, the r A S scheme is dealing with particular 
solutions less rapidly decreasing. 
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FIG. 6. Error on the resolution of the scalar Poisson equation for a solution containing 
bounded logarithm functions. 
The scale for the number of coefficients is linear. 

The circles represent the error in the kernel, the squares in the shell and the diamonds in the 
external domain. 
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FIG. 7. Error on the resolution of the scalar Poisson equation for a solution containing In 
functions for I = 2. 

The scale for the number of coefficients is logarithmic. 

The circles represent the error in the kernel, the squares in the shell and the diamonds in the 
external domain. 

Solid lines represent the scheme with r 4 S, the dotted ones the scheme with r 3 S and dashed lines 
the one with r 2 S. 
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FIG. 8. Exponent of the power-law followed by the error shown in Fig. 7 as a function of 
the index /. 

The solid lines correspond to the r 4 S scheme, the dotted lines to the r 3 S one and the dashed lines 
to the r 2 S one. 

The circles represent the error in the kernel, the squares in the shell and the diamonds in the 
external domain. 



In Fig. 8 the slope of the power-law is plotted as a function of the harmonic 
index I, for the three different schemes. It reveals an error decreasing as 7V~ 2 (' +1 ) 
for the r 2 S and r 3 S schemes and as N~ 21 for the r 4 S one. Let us mention that the 
r 2 S scheme yields an error following the same power-law than the one rising from 
the description of the associated function (cf. Sec. 3.1.3) , making us confident 
about the origin of such a behavior. 



3.3. Convergence properties 

All the examples shown in the previous section enable us to propose the two 
following empirical properties concerning the decrease of the error. 

• Property 1 : if the source is decreasing as r~ k at infinity and does not contain 
any spherical harmonics with I > k — 3, then the error is evanescent. 

• Property 2 : if the source decreases at least as r~ k at infinity, then the error 
decreases at least as N~ 2 ( k ~ 2 ^ (rcsp. N~ 2k ) for the r 2 S and r 3 S schemes (resp r 4 S 
scheme) . 

The first one is just issued from the condition to have a In function in the external 
domain and the second one from the values of the power-law found in the previous 
section. 



20 



P. GRANDCLEMENT ET AL. 



4. VECTORIAL POISSON EQUATION 

Using the Poisson equation solver from Sec. 2 and studied in Sec. 3, we focus 
now on the vectorial Poisson equation given by Eq. (2), in the non-degenerated 
case (i.e. A ^ —1). 

Let us first mention that the operator A + AV(V- ) has been shown to be strongly 
elliptic and self-adjoint in Refs. [5, 6] in the case A = 1/3 (conformal Laplace 
operator). Conditions for existence and uniqueness of solutions have been presented 
in Appendix B of Ref. [4] . The harmonic vectorial functions of this operator and 
the associated multipole expansions have been discussed by O Murchadha [21]. 

Three different schemes have been previously proposed by other authors [14, 15, 
16] to reduce the resolution of Eq. (2) to those of four scalar Poisson equations. Let 
us emphasize that those three schemes are not covariant. They are only applicable 
in Cartesian coordinates which allow us to commute operators like Laplacian and 
gradient. 

Let us mention the fact that a different method, based on solving for the degen- 
erated case (i.e. A = 1) has been proposed in [10] but is not studied in the present 
work. 



4.1. The Bowen-York method 

The idea of this method (sec [14]) is to search for the solution of Eq. (2) in the 
form 

N = W + V X , (44) 

where W and \ are solutions of 

AW = S (45) 
A X = --^-yV- (46) 



This method gives a solution to Eq. (2) but let us check that this solution is the 
one that is C . W is C , being solution of a Poisson equation. This implies that 
the source of the equation for x is continuous, and that x is C 2 . This is sufficient 
to ensure that N is C 1 . The scheme finds the only solution C 1 and going to zero at 
infinity. 

Unfortunately this very simple method is not applicable with our Poisson equa- 
tion solver, for the physical sources are not decreasing fast enough at infinity. For 
the problem that motivated this study, namely binary neutron star systems [22, 23], 
the source S of Eq. (2), is expected to behave like r~ 4 at infinity implying that 
we can calculate W. This vector field is acting like r _1 at infinity, for r _1 is an 
homogeneous solution of the scalar Poisson equation usually present (monopolar 
term). 

So the source of the equation for x, being the divergence of W, behaves like 
r~ 2 . This decreasing is not fast enough to compute the value of x- Analytically 
no problem occurs because only the gradient of x is relevant, not \ itself, for 
the calculation of the solution. To summarize, the implementation of this scheme 
conducts to the computation of diverging quantities, making the result wrong in 
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the external domain. We should say that this scheme is applicable for domains not 
extending to infinity. However it may be possible to use it by treating analytically 
the diverging quantities. 

4.2. The Oohara-Nakamura method 

In this case (see Sec. 3.1.1 of [15]) we start by solving the following scalar equation 

A X = • S. (47) 

Then the solution of Eq. (2) is found by solving the following set of three equations 

AN = S-\V X - (48) 

Comparing (2) with (48) shows that this scheme gives the exact solution of Eq. 
(2) if and only if 

Vx = V (V • ffj . (49) 
But the scalar equation (47) only ensures that 

A (x - V • Nj = 0. (50) 

From the study of Sec. 2, wc can show that it is possible to construct an homo- 
geneous solution of the scalar Poisson equation, in all space, that is non-zero, going 
to zero at infinity, if and only if, that solution is not C . 

In the general case, V • N is only C° at boundary between the different domains, 
while x, solution of a Poisson equation, is C x . So it is possible to fulfill Eq. (50) 
and not Eq. (49). If V • N is C 1 , then Eq. (50) implies, as shown by Sec. 2, that : 
X = V • N. In this case, the condition (49) is trivially fulfilled. Imposing that V • N 
is at least C 1 , is equivalent to impose that S is continuous across every boundary. 

To conclude, let us say that the Oohara-Nakamura method gives the exact so- 
lution if and only if the source S is continuous across every boundary delimiting 
the different domains. This property is general, meaning that it is not due to our 
numerical method. We can mention that the found solution is the C 1 one, because 
it is calculated as solution of three scalar Poisson equations. 

Second let us see if this scheme is applicable, using our scalar Poisson equation 
solver. At first sight, this scheme suffers the same drawback than the Bowen-York 
one. Because of homogeneous solutions of scalar Poisson equation, x is decreasing 
as r _1 at infinity and its gradient as r -2 , which is not enough to allow us to solve 
the set (48) of three scalar Poisson equations. 

The difference is that the solution of Eq. (48) is the solution of the vectorial 
Poisson equation (2) and we must be able to set it to zero at infinity, contrary to 
the Bowen-York method where the problem occurs for auxiliary quantities. 

So it must be possible to show that the source of Eq. (48) decreases fast enough, 
that is, at least as r -3 . The problem arises from the monopolar term of \i i- e - t ne 
only one that gives an homogeneous solution in in the external domain. It is 
known, that the mono-polar term M of the solution of a scalar Poisson equation 
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with source a, is given by 

Mo = h SIS ffd3r ' (51) 

the integration being performed over all space. 

Now we have a = V • S. The use of Green formula leads to 

M o = iff S-ds, (52) 

the surface integration being done at infinity. But S decreases as r~ 4 , implying 
that the surface integral is zero, that is M = 0. This remains true if the source 
acts only like r~ 3 . 

So the monopolar term of \ is zero, which implies that X decreases at least as 
r~ 2 . This behavior ensures that the source of Eq. (48) decreases as r -3 , allowing 
us to find the unique solution going to zero at infinity. 

We implemented and tested this scheme, recalling the reader that it is only ap- 
plicable if the source of Eq. (2) is continuous and requires that the source decreases 
at least like r~ 3 at infinity. 

4.3. The Shibata method 

The solution is now searched as (see [16]) 

N=l^ 1 W- 1 - J ^- 1 (v X + VW.f), (53) 
where W and X are solutions of 

AW = S (54) 

A X = -r-S. (55) 

and r denotes the vector of coordinates (x, y, z). 

This scheme gives a solution to Eq. (2), but, as for the Bowen-York method, let 
us quickly check that it is the unique C 1 going to zero at infinity. At infinity, W, 
solution of scalar Poisson equation, is behaving at least like r _1 . This ensures that 
VVT • r is zero at infinity, proving that the solution goes to zero. 

Concerning the continuity, being solutions of scalar Poisson equations, we know 
that both W and X are at least C x . But we have to take care about the term 
+ VW ■ r of Eq. (53). First we can show that 



A(r-W) = f-S + 2V-W. (56) 

Using that property and the equation for X we can see that 

A{r-W + x) =2V-W. (57) 

The source of that equation is C°, so that r ■ W + X is C 2 . The term of Eq. (53), 
can be expressed as 

Vx + VW>-f= v(r-W + x) -W. (58) 
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Using the continuity properties found above, it is easy to see that the right-hand 
side of Eq. (58) is C 1 , which ends our demonstration by proving that the calculated 
N is C 1 . 

As before, let us now check if this method is applicable by means of our scalar 
Poisson equation solver. The source of the equation for \ decreases at least like r -3 
at infinity if and only if S decreases like r~ A . Like the Oohara-Nakamura scheme, 
this one does not involve any diverging quantities and so is suitable for numerical 
purposes. 

This method has been implemented and, contrary to the Oohara-Nakamura 
method, can be used even with discontinuous source, but requires that S decreases 
at least like r -4 at infinity, which, let us recall, is the case for the physical problems 
we intend to study. 

4.4. Convergence criterion 

As seen before the resolution of Eq. (2) reduces to that of four scalar Poisson 
equations. So we should be able to use the results of Sec. 3.3 to established 
convergence criterion for the schemes proposed in [15, 16]. 

4-4-1. The Oohara-Nakamura scheme 

Let us suppose that the source S of Eq. (2) contains only one spherical harmonic 
Y" ; m and decreases as r~ k at infinity (k > 3). 

For the Oohara-Nakamura method, the source of the first Poisson equation is 
V • S : the degree of the harmonic is I + 1 and the decrease is as r~( fe+1 ) . These two 
effects are opposed concerning the convergence properties established in Sec. 3.3. 
So, in the case no logarithm appear during the calculation to find x, x contains one 
spherical harmonic I + 1 and decreases as r~( fe_1 ) and so Vx, part of the source of 
Eq. (48), contains one spherical harmonic with I + 2 and acting like r~ k at infinity. 
So the conditions for the appearance of a Gibbs-like phenomenon are "harder" 
by two degrees than for a scalar Poisson equation and occurs for a source with a 
spherical index 1 + 2. 

4-4.2. The Shibata scheme 

Suppose we consider the same source than in the previous section. The conver- 
gence properties for the equation for W are the same than for a usual scalar Poisson 
equation. 

Concerning the equation for \ the source is — r ■ S. Performing such an operation 
on S increases the degree of the spherical harmonics by one unit. At the same time, 
the decrease of the source is slower, due to multiplication by r everywhere. Those 
two phenomena have the same effect on the convergence criterion we previously 
established. As for the Oohara-Nakamura scheme, the criterion are "harder" by 
two degrees but the Gibbs-phenomenon occurs for a source in I + 1 . 

4-4-3. Convergence properties 

We are now able to deduce convergence properties for the two schemes. From the 
study above, we can see that if the condition for the appearance of the Gibbs-like 
phenomenon is the same, it is not associated with the same index /. This results in 
the two following properties : 



24 



P. GRANDCLEMENT ET AL. 



• Property 1 : if the source of a vectorial Poisson equation is decreasing as 
r~ k at infinity ( k > 3 for the Oohara-Nakamura scheme and k > 4 for the Shibata 
one) and does not contain any spherical harmonics with I > k — 5, then the error 
is evanescent. 

• Property 2 : if the source decreases at least as r~ k at infinity then the error 
is decreasing at least as _/V~ 2 ( fc ~ 2 ) for the Oohara-Nakamura method (k > 3) and 
at least as N~ 2( - k ~ 3 ^ for the Shibata one (k > 4). 



5. ACCURACY OF THE VECTORIAL POISSON EQUATION 
SOLVERS ESTIMATED BY COMPARISON WITH ANALYTICAL 

SOLUTIONS 

To check the validity of the schemes and their convergence, we used the same 
method than for the scalar Poisson equation, that is the use of analytical solutions 
of various properties. The solutions associated with the sources have been obtained 
by following analytically the Shibata scheme. 

5.1. Continuous source 

Let us consider the case of a continuous source extending to infinity, say for 
example, in the external compactified domain, for r > R 

and for r < R 

S * = ' SV = R^ ' S * = RP+S' ^ 

Note that this source is C°, the minimum requirement for the Oohara-Nakamura 
method to be applicable. 

For fi^O, the associated solution in the external domain is 

N x = 1 X n + 5 x (61] 

(A + l)n(n + 3)r"+ 3 (A + 1) 15n R n r 3 { ' 



and for r < R 



9 

XT 



N x = ' (62) 

10(A + l)i?»+ 5 (A + 1) 6 (n + 3) R n + 3 K ' 

the other components being obtained by permutation of x, y and z. 

For n / no Gibbs-phcnomenon occurs by solving the equations with S as 
source. For n < 2, a Gibbs-like phenomenon should appear due to the vectorial 
nature of Eq. (2). It is not the case because of simplifications due to the symmetry 
of the source. It just shows that the two convergence criteria established above are 
rather pessimistic. The evanescent error is shown in Fig. 9. As for the scalar case, 
a saturation is attained at a level of approximatively 10 -11 . 

5.2. A vectorial Gibbs-like phenomenon 

At this point, we would like to exhibit an analytical solution that produces a 
Gibbs-phenomenon which arises from the vectorial nature of Eq. (2). Let us con- 
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FIG. 9. Error on the z component for a continuous source extending to infinity (Eqs. 59-60 
with n = 1). 

The scale for the number of coefficients is linear. 

The solid lines represent the Shibata scheme and the dashed lines the Oohara-Nakamura one. 
The circles represent the error in the kernel, the squares in the shell and the diamonds in the 
external domain. 



sider the following source 



S z = 



(63) 



in the external compactified domain, and for r < R 



S z = 



R 7 ' 



(64) 



We set the two other components to zero in all space. 

If we solve the scalar Poisson equation with S z as source, the error will be evanes- 
cent, as shown by the study of Sec. 3.3. But, according to the conclusion we 
obtained concerning the convergence criterion of a vectorial Poisson equation, a 
Gibbs-like phenomenon should appear due to the vectorial nature of Eq. (2). 
In the external domain, the associated solution is 



N y 



1 A 



2 A 



r 2 

Z X 



- 1 

59 

r u \350 
1 A 
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V 
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350 



+ 



14 
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(65) 
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FIG. 10. Error on the z component for a source implying a Gibbs-like phenomenon. 
The scale for the number of coefficients is logarithmic. 

The solid lines represent the Shibata scheme and the dashed lines the Oohara-Nakamura one. 
The circles represent the error in the kernel, the squares in the shell and the diamonds in the 
external domain. 
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and for r < R, we found 



N x 



N z 
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1 A 
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z 
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As expected, Fig. 10 shows an error obeying a power-law. This feature is more 
evident in the external domain where the particular solution is directly present. 
The Gibbs-like phenomenon appears for the two schemes. Let us apply Property 2 
to determine the exponent of the power law. The source of the equation decreases as 
r~ 6 . This implies that the error for the Oohara-Nakamura scheme should decrease 
at least as N~ 8 and as ./V -6 for the Shibata one. This is well confirmed for the 
Shibata scheme which exhibits an exponent —6.4. For the Oohara scheme it turns 
out that the criterium is rather pessimistic for the error decreases faster than A^ -12 . 
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5.3. A discontinuous source 

As previously explained, the Oohara-Nakamura scheme fails to solve Eq. (2) in 
the case of a discontinuous source. We will now consider such a source and show 
that the Shibata method is efficient, even in such a case. 
In the compactified domain, r > R, we choose the following solution 

N x = — (67) 

For r < R, we ensure the continuity of the solution and its derivative by choosing : 

N x — x (ar 6 + br 4 ) , (68) 

A + n 6 + n m1 . . 

where a = — and b = — . I he associated source is obtained by calcu- 

lating the left-hand side of Eq. (2). In the external domain we obtain 

S- = n(n-3-3X)-^ + n(n + 2)X-^ (69) 

2 

S y = -Xn^—+n{n + 2)X^j 

-n+2 v ' r n+4 



r 

S z = -\n—— +n(n + 2)X 



Z • , x 2 z 



and for r < R, we have 



S x = x [(54 + 18A) ar 4 + (28 + 12A) br 2 ] + Xx 3 (24ar 2 + 8b) (70) 
S y = \y [6ar A + 4br 2 + x 2 (24ar 2 + 8b)] 
S z = Xz [6ar 4 + Abr 2 + x 2 (2Aar 2 + 8b)] . 



Depending on the value of n, the error is evanescent or not. Only a few spherical 
harmonics are present in the source and we can show that we expect, for example, 
an evanescent error for n = 4 and a Gibbs phenomenon for n = 5. This might seem 
not to be in agreement with the convergence criterion previously established, but 
the reader should recall that they are rather general and much more pessimistic to 
handle simple sources such as the ones which are considered here. 

The results presented by Figs. 11 and 12 show that the discontinuity of the 
source has no effect on the resolution of the vectorial Poisson equation, as long as 
the Shibata scheme is used. For n = 5, the source is like r~ 6 at infinity and we 
expect an error decreasing more rapidly than iV~ 6 . Fig. 12 shows an extremely 
good agreement with the prediction, for the power-law exhibits an exponent of 
-6.4. 

6. DEVELOPMENTS 

In this section we present some extension of this work, that is solving more 
complicated equations using the schemes presented here as milestones. 

The first extension that has been conducted regards non-spherical domains, with 
spheroidal shapes (i.e. they must have the same topology as a sphere). This is 
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FIG. 11. 

Error on the x component for discontinuous source (Eqs. 69- 70 with n = 4). 
The scale for the number of coefficients is linear. 

The circles represent the error in the kernel, the squares in the shell and the diamonds in the 
external domain. 
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FIG. 12. Same as Fig. 11 but for n = 5 ; the scale for the number of coefficients is now 
logarithmic. 
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very useful for we can define the boundary of each domain to match with surfaces 
of discontinuity, like stellar surfaces, so that each field is C°° in each domain pre- 
venting any Gibbs phenomenon. Thanks to some mapping onto a sphere, solving 
the Poisson equations with such boundaries reduces to the spherical case, with cor- 
rection terms appearing in the source. The equation is then solved by iteration. 
The method is described in details in [13]. In that paper the calculation of the 
Mac-Laurin and of the Roche ellipsoids have been compared with the analytical 
solutions. The behavior of the error when one increases the number of coefficients 
happens to be evanescent (see Figs. 5 and 6 of Rcf.[13]). Those calculations being 
done in the Newtonian case, all the sources are compactly supported. This shows 
that the non-sphericity does not introduce any new Gibbs phenomenon with respect 
to the spherical case. 

Concerning calculations in general relativity (i.e. with sources extending to in- 
finity), results have been obtained for rapidly rotating strange stars in [20] using 
non-spherical domains. Convergence properties have not been fully explored for 
there exist no analytical solution to compare with. Anyway we can suppose that, 
the sources containing almost every spherical harmonics, the convergence will no 
longer be evanescent but rather follows a power-law. 

Another important extension of this work is to deal with two bodies, for example 
orbiting binary neutron stars. This case has been successfully studied in [22, 23], 
by means of the Poisson solvers presented here. The main difference with the cases 
we discussed in the present paper is that the sources are no longer spheroidal but 
are concentrated on two spheroidal domains, being the two stars. An equation of 
the type (1) is then split into two parts : 



where the real source is S = Si + S 2 - We use two sets of spherical coordinates, 
one centered on each star and the splitting is done so that Si is mainly centered 
on the first star and S2 on the other one (see [23] for details). The sources Si are 
then well described in spheroidal topology and the total equation is well solved, 
the solution being F = F\ + F 2 . We used that method to compute Newtonian 
configurations and compare them with semi-analytical solutions. Fig. 13 shows the 
error done for the same configuration as Fig. 7 of [23] for a coordinate separation 
of 100 km. This calculation being Newtonian, the sources are compactly supported 
and the error seems to be evanescent, but we have to be cautious for the number 
of coefficients of the expansion arc not maintained fixed. Extensive convergence 
properties have not been conducted but it seems that the splitting of the equation 
into two parts does not introduce any new Gibbs phenomenon. As for the single 
body problem, convergence of calculations with sources extending to infinity (i.e. 
in general relativity), have not been studied but we expect a Gibbs phenomenon to 
occur for the sources contain almost every spherical harmonics. 

To finish with the extension of this work let us mention the case of black holes. In 
that case the equations are not solve in all space but only on the domain exterior 
to the holes horizons. This means that we have to remove the kernel from the 
computational domain. The regularity condition at the origin is then replaced by 



AF 2 



- Si 
= S 2 



(71) 
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FIG. 13. Relative error, estimated by means of the virial theorem, for a Newtonian 

irrotational binary star calculation with respect to the number of Chebyshev coefficients. 

a boundary condition on the boundary of the innermost shell. We have been able 
to use that to impose condition on the value of the solution (Dirichlet problem) or 
on its first radial derivative (Neumann problem) . This extension has nothing to do 
with the compactified domain and we expect the convergence properties to be the 
same as the one exhibited in the present work. We are currently applying this to 
compute realistic physical binary black holes configurations. 

7. CONCLUSION 

We have presented a scalar Poisson equation solver based on spectral method. It 
enables us to solve the Poisson equation for a source extending to infinity and going 
to zero at least like r~ 3 . Our multi-domain approach enables to deal with a source 
which is C°° in each domain. Nevertheless some Gibbs phenomenon can appear 
due to the existence of particular solutions which contain logarithm functions in 
the external domain. Such functions are not well described in terms of Chebyshev 
polynomials, resulting in a Gibbs-likc phenomenon. We exhibited the conditions 
for the appearance of such an effect and quantified it, leading to the conclusion 
that, for a source decaying as r~ k (k > 3), the error of the numerical solution 
is evanescent if the source does not contain any spherical harmonics with index 
I > k — 3. Otherwise, the error decreases at least as N~ 2 ( k ~ 2 ^ , N being the number 
of Chebyshev coefficients. 

We used this scalar Poisson equation solver to solve the generalized vectorial 
Poisson equation given by Eq. (2) for a source going to zero at least like r~ A . 
Three different schemes have been discussed. We showed than the one proposed by 
Bowen & York [14] is not applicable to domains extending up to infinity, by means 
of our methods, because it gives rise to diverging auxiliary quantities. The scheme 
from Oohara & Nakamura [15] is applicable as long as the source is continuous and 
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has been successfully implemented. The last scheme, proposed by Oohara, Naka- 
mura and Shibata [16], is applicable even for discontinuous sources and has been 
successfully implemented too. The convergence properties of the two implemented 
schemes have been derived from the ones of the scalar Poisson equation solver and 
checked by comparison between calculated and analytical solutions. 
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